module FullLinearSolver

    real(kind=8), parameter :: errorLimit = 1.0d-14

contains

    function GaussianEliminationAlgorithm(A) result (r)
        implicit none

        real(kind=8), dimension(:,:), intent(inout) :: A

        real(kind=8), dimension(:), allocatable :: r

        integer :: n
        integer :: k, i, j

        real(kind=8) :: base
        n = size(A,1)

        if (n == 0) return

        allocate(r(n), stat=k)

        if( k /= 0) then
            print *, 'Error allocating memory for response vector.'
            return
        end if

        do k = 1, n-1
            if ( abs(A(k,k)) <= errorLimit ) then
                print *, 'The Matrix presents stability problems. Exiting Gaussian Elimination.'
                return
            end if
            do j = n+1, k, -1
                A(k,j) = A(k,j)/A(k,k)
            end do
            do i= k+1, n
                if( A(i,k) == 0 ) continue
                base = - A(i,k)
                do j = k, n + 1
                    A(i,j) = A(i,j) + base * A(k,j)
                end do
            end do
        end do
        r(n) = A(n, n+1)/A(n,n)
        do k = n - 1, 1, -1
            r(k) = A(k, n+1)
            do j = k + 1, n
                r(k) = r(k) - A(k,j) * r(j)
            end do
        end do
    end function


    subroutine ShowMatrix(A)
        implicit none

        real(kind=8), dimension(:,:) :: A

        integer :: m, n
        integer :: i, j

        m = size(A,1)
        n = size(A,2)

        do i = 1, m
            do j = 1, n
                write(*,'(F13.5)',advance="no") A(i,j)
            end do
            write(*,*)
        end do

    end subroutine

    subroutine LinePivoting(A)
        implicit none
        real(kind=8), dimension(:,:), intent(inout) :: A

        integer :: m, n
        integer :: k, i, j
        real(kind=8) ::swap

        m = size(A,1)
        n = size(A,2)
        do k = 1, m-1
            do i = k+1, m
                if ( abs(A(k,k)) < abs(A(i,i)) ) then
                    do j= 1, n
                        swap = A(k,j)
                        A(k,j) = A(i,j)
                        A(i,j) = swap
                    end do
                end if
            end do
        end do

    end subroutine

end module
